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We present a tutorial discussion of Monte Carlo methods for equilibrium and nonequilibrium 
problems in interfacial electrochemistry. The discussion is illustrated with results from simulations 
of three specific systems: bromine adsorption on silver (100), underpotential deposition of copper 
on gold (111), and electrodeposition of urea on platinum (100). 

I. INTRODUCTION 

In the context of interfacial electrochemistry the bulk electrolyte can conveniently be thought of as a large, homo- 
geneous reservoir of charge, adsorbate particles, and energy. Then one can focus on the metal-electrolyte interface, 
which can be viewed as an essentially two-dimensional system of ions and molecules at or near a solid surface. Conse- 
quently, interfacial electrochemistry can profitably be thought of as a branch of surface science. It is therefore natural 
that a number of experimental and theoretical methods originally associated with more traditional branches of surface 
science now are routinely applied in interfacial electrochemistry, including several that are discussed in this volume. 

On the theoretical side, one such method is computational lattice-gas modeling. This method is particularly well 
suited to the study of fluctuations and phase transitions in systems characterized by specific chemisorption at well 
defined sites on the surface. Different aspects of this approach have been reviewed several times by others as well 

as by one of the present authors 1^-0]. Here, it will only be described briefly. Monte Carlo simulation of lattice-gas 
models as a tool for understanding interfacial electrochemistry is the main topic of this Chapter. 

Why bother with Monte Carlo simulation when there are so many mean-field techniques, ranging from the simple 
Frumkin and Temkin isotherms to sophisticated cluster-variation methods? This is a natural question, and there are 
three compelling answers: 

1. Mean- field methods are really not very accurate for most two-dimensional systems 1^. This is not to say that 
very useful results cannot be obtained by mean-field methods, for they certainly can, but simple mean-field 
methods ignore the spatial structure of the system. A nontrivial amount of work and attention goes into 
manually incorporating enough of this structure to give accurate results, especially in the vicinity of phase 
transitions. In contrast, Monte Carlo codes naturally take spatial fiuctuations into account. Most of these codes 
are quite simple to write and easily adapted to model different systems. 

2. Monte Carlo methods are not limited to the study of equilibrium properties. Dynamic Monte Carlo simulations 
provide a tool to study time dependent phenomena in systems where spatial fluctuations may strongly modify the 
time evolution. Applied in this way, Monte Carlo simulation provides a more detailed supplement to standard 
rate-equation methods. 

3. Computer technology has now evolved to a point where extensive simulations can be performed on work stations 
or personal computers that are well within the budget for most groups. Actually, several computers that 
are perfectly able to handle large Monte Carlo simulations may already be sitting in your laboratory! The 
competitive edge of methods that are analytically solvable or rely on numerical solution of relatively small sets 
of equations is much smaller than it was only a few years ago. 

The purpose of this Chapter is to discuss some of the basic aspects of Monte Carlo simulation directly in the context 
of interfacial electrochemistry. It is not intended to be a full textbook on Monte Carlo simulation. There is no sense 
in duplicating the excellent textbooks PJTc|] and large tutorial articles ||l^,|l^ that already exist. Rather, we hope to 
illustrate some of the details and pitfalls associated with constructing and running Monte Carlo simulations of specific 
systems of interest in surface electrochemistry. 

The outline of the rest of this Chapter is as follows. The lattice-gas approach to modeling electrochemical adsorption 



is briefly introduced in Sec. [I A. The process of comparing numerical and experimental results to adjust the lattice- 



gas parameters is the subject of Sec. II B, The basics of equilibrium and dynamic Monte Carlo simulations as 



nonperturbative methods for analyzing lattice-gas models at finite temperature are discussed in Sec. [II. These 
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simulations can be used to obtain numerical results for such experimentally measurable quantities as adsorption 
isotherms, voltammetric currents and charge densities, and images obtained by either diffraction techniques or atomic- 
resolution microscopies. In Sec. 0, we discuss three illustrative examples drawn from our own research experience: 
the adsorption of bromine onto Ag(lOO) ||l^-[l5|], the underpotential deposition (UPD) of copper on Au(lll) from 
a sulfate-containing electrolyte P,p^-pl| , and the electrosorption of urea on Pt(lOO) from an acid electrolyte [|6|,p2[. 
Particular emphasis is given to the simulational aspects of these projects. Since this is not a review article, references 
to the original literature have been limited to those we thought most illustrative for our purposes. 



II. LATTICE-GAS MODELS OF CHEMISORPTION 



A. Structure of the Models 



In lattice-gas models adsorption is only allowed on a discrete lattice of adsorption sites. The concentration of 
adsorbate at site i, Ci, is 1 if the site is occupied by the adsorbate and otherwise. The energy of a particular 
configuration of occupied sites, c, is given by the grand-canonical lattice-gas Hamiltonian 



n{c) = -Ai E ^» - E E ^^^^ - ^3(c) . (1) 




The electrochemical potential of the adsorbate in solution is /l, and the first term on the right-hand side is the free 
energy of adsorption when adsorbate-adsorbate interactions are ignored. The next term represents all of the two-body 
adsorbate-adsorbate interactions. The index n indicates the separation between sites, e.g. n = \ indicates nearest 
neighbors, while indicates the sum over all pairs of neighbors of rank n. The interaction energy for two adsorbate 

particles occupying n-th nearest neighbor sites is —$("). Interactions between three or more adsorbed particles are 
possible [^,^ — here they are symbolically represented as the term Tis. The sign convention is such that /2>0 denotes 
a tendency for adsorption in the absence of lateral interactions, and $^"^>0 denotes an effective attraction. Models 
of particular electrochemical systems differ in their binding-site geometries and the values and ranges of the lateral 
interactions. 

The connection of the theoretical lattice-gas model to the chemical and electrochemical variables that characterize 
the experimental system is made through the electrochemical potential. The electrochemical potential of X is related 
to its bulk activity [X] and the electrode potential E as 

-^ = ^° + RT\n^-zFE , (2) 

where R is the molar gas constant, T is the absolute temperature, F is Faraday's constant, and z is the effective 
electrovalence of X. The quantities superscripted with a are reference values which contain the local binding energies. 
Here we treat them as adjustable model parameters. 

The thermodynamic density conjugate to the electrochemical potential is the surface coverage 



e = 7v-i^c, , (3) 



where N is the total number of adsorption sites. With the sign convention that oxidation currents are positive, the 
adsorbed charge per adsorption site is 

q = -ezQ , (4) 

where e is the elementary charge unit. 

Many interesting systems involve the cooperative chemisorption of two different species, say A and B, onto the 
surface. In this case there is a concentration for each species at site i, and cf , with simultaneous adsorption 
of both species at one site forbidden . The lattice-gas Hamiltonian for such a two-component system can be 

written as |^ 



U (c) = ~^1A E ~ E ' 

^3 . (5) 



E 



(n) (n) (n) 

^'aI E + '^'bI E ^"s" + E (^^" + ^^4) 

(ij) (ij) (ij) 
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The total configuration of the system, which depends on both the configurations of A and B particles, and c^, is 
still represented as c. The interactions ^xy (where X and Y can be A or B) now depend on the adsorbed species as 
well as the distance between sites. The higher-order interactions, Ti.3, may also involve combinations of A and B. 

Each component has its own electrochemical potential fix, surface coverage 0x, and effective electrovalence zx- 
When the temperature is fixed, a phase diagram can be constructed in the (/IajA'b) plane. For a specific chemical 
system with fixed bulk activities [A] and [B], Eq. (H) for /Ia and /2b form the parametric representation of a line in 
the phase diagram. The line has slope za/^b, and the electrode potential E is the parameter that determines the 
location along this isotherm. 

Since the entropy vanishes at zero temperature, the T = phase diagram, or ground-state diagram, can be calculated 
analytically using Eq. (^. In principle this can be done directly, but in practice another approach is often used. For 
a given point in the phase diagram, the ground state corresponds to the phase with the lowest energy. A list of 
the conceivable ground states compatible with the lattice symmetry [|6| and their energies is generated. Then a 
ground-state diagram can be constructed automatically by finding the lowest-energy phase directly while scanning in 
the p,A and /2b directions. 

In contrast, the finite-temperature phase diagram for most systems can only be mapped out using numerical 
methods, such as the Monte Carlo simulations that are the main subject of this Chapter. Since the number of 
possible ground states for systems with longer-range interactions can be very large, it is in practice advisable to check 
the ground-state calculations and Monte Carlo simulations at low temperature against each other for consistency. 

When two types of ions adsorb onto the surface, the adsorbed charge density is 



q ■ 



-e (2;a6a + ^^bBb) 



(6) 



In the absence of diffusion and in the limit of vanishing potential sweep rate , the voltammetric current can be 
estimated in terms of the equilibrium lattice-gas response functions dQx/dp-Y- Specifically, using the chain rule this 
quasi-equilibrium current is obtained as 



dq 
dt 



dq dE 
dE~dt 



The total differential of the charge density is 

dQA 



dg = 



ZA 



dfjLp 



ZB 



dfiA + ZA 



ZB 



6>9e 

dfj-B 



dfiB 



(7) 



(8) 



which can be simplified with the Maxwell relation 90a/9/2b — 90b/c^/2a- For aiw particular electrochemical ex- 
periment, d/2x = —zxFdE. Substituting this into Eq. (j^ and that result into Eq. (^, the quasi-equilibrium current 
becomes 



eF 



za 



OQa 
dfiA 



2zazb 
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+ zt 



dfjLB 



dE 
dF 



(9) 



The response function can, in turn, be estimated from the equilibrium partition function using standard statistical- 
mechanics arguments, as described in the next section. 

Coverages, charges, and currents are given here per adsorption site. They can be converted to the experimentally 
measured units of inverse surface area through division by the area per adsorption site. 



B. Model Development and Refinement 



Accurate theoretical description of lateral interactions in an adsorbed layer must consider the complete substrate- 
adlayer- solution system. For instance, the geometry of adlayers observed in situ for CO on Pt(lll) are different for 
ultra- high vacuum, high-pressure gas, and electrochemical environments p8| . In addition, adsorption or solvation at 
one point on the surface can affect adsorption at another point through deformation of the substrate or distortion of the 
substrate electron structure. Because of the complicated nature of these systems, deriving accurate ah initio estimates 
of interaction energies in electrochemical adsorption problems is presently not feasible. A practical alternative is to 
treat the interactions as adjustable parameters. Because of the large number of parameters and the time-consuming 
simulations, parameter estimation by a formal optimization procedure must often be abandoned in favor of adjusting 
the parameters manually to incorporate experimental results and other chemical information efficiently. This method 
also has its difficulties. In particular, all of the potential consequences of adjusting a parameter may not be easily 
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foreseen, and often several parameters must be adjusted in concert to change a single aspect of the model. In addition, 
there is no guarantee that the observations from two radically different sets of interactions are measurably different 
in the region explored by the simulations. Nevertheless, the encouraging results of previous lattice-gas studies of 
electrochemical systems |6|,^- ^^5|j29|j30[| indicate that when proper attention is paid to including all available 
experimental information consistently, this approach has considerable predictive power. Furthermore, as effective 
interactions obtained by first-principles calculations become available in the future, results from lattice-gas models 
will provide crucial information for testing the consistency of such first-principles interactions with experimentally 
observed thermodynamic and structural information. 

The modeling strategy starts by using theoretical and experimental knowledge about the substrate structure and 
the shapes and sizes of the adsorbate particles to formulate a specific lattice-gas model. The coverages and adlayer 
structure of different adsorbate phases can be determined from microscopy, scattering, radiochemical experiments, or 
standard electrochemical methods. Group-theoretical ground-state calculations |2^] are then performed to determine 
a minimal set of effective interactions consistent with the observed adsorbate phases A ground-state 

diagram is obtained by pairwise equating the ground-state energies of the different phases, as obtained fro m the 
lattice-gas Hamiltonian, e.g. Eq. (|l]) or Eq. (^). A typical example, corresponding to the model discussed in Sec. [VB, 
is shown in Fig. ^. 

The finite-temperature properties of the lattice gas, calculated by Monte Carlo or other methods, [^-| 19 2^j25| , p^ , ^ 
can be compared to available experiments and used to refine the effective interactions. The process should be iterated 
until satisfactory agreement between model and experiment is achieved. New experiments may be required to clarify 
the properties of the chemical system and refine the theoretical model. 



III. MONTE CARLO SIMULATION 

A. Monte Carlo Sampling of Equilibrium Configurations 

In statistical mechanics the properties of a system in equilibrium, such as the coverages and response functions 
mentioned in the previous section, are calculated from the partition function 

Z = ^exp[-H(c)/i?T] , (10) 

where J2{c} '^^^^ possible configurations of the adsorbates. The expectation value of a quantity Y at equilibrium 
is the weighted sum 

{Y)^jY.Y{c)exp[~'Hic)/RT] , (11) 

{c} 

where the weight is a Boltzmann factor. For most realistic models it is hard to calculate, or even estimate, Eq. (|l0| ) 
or (|l^) analytically. On the other hand, Monte Carlo simulation methods can be used to estimate the properties of 
even quite complicated models numerically 1^. 

The reason for our advocacy of the nonperturbative numerical Monte Carlo approach is that such calculations are 
much more accurate for two-dimensional systems than even quite sophisticated mean- field approximations From 



12 1, using modest computational 
121, Monte Carlo methods have 



Monte Carlo simulations one can obtain thermodynamic and structural information | 
resources and programming effort. In combination with finite-size scaling analysis 
contributed significantly to the theoretical understanding of fluctuations and ordering in adsorbate systems. Monte 
Carlo codes also have the advantage that they are relatively easy to modify to accommodate changes in lattice 
structure and/or interaction geometries and ranges. This is useful when one wants to modify the model to improve 
the agreement with experiments, or refurbish an existing code to model a different system. 

One Monte Carlo method for estimating Eq. (^l]) is to evaluate the sum directly from randomly generated adsorbate 
configurations. For most problems where temperature is important, this method is inefficient since only the most 
energetically favorable configurations make significant contributions to the sum. A much more efficient method, called 
importance sampling, samples more frequently from the energetically favorable configurations. Usually this is done 
by using the Boltzmann weight, the probability of an equilibrium system being in a particular state, to determine the 
sampling probability of that state. In addition, one can use the idea of detailed balance, that the rate of going from 
state c to state c' in equilibrium is the same as the rate of going in the opposite direction. If M{c) is the density of 
states for configuration c and the probability of going from configuration c to configuration c' during a certain period 
of time is TZ{c c'), the transition rate is N{c)TZ{c c'). Detailed balance implies that 
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AA(c) ^ n{c' ^ e) 

N{d) n{c ^ c') ■ ^ ^ 

To sample configurations according to tlie Boltzmann weights, the transition probabihties need to obey 

^(^'-^)-expr-M2)^V (13) 



An efficient Monte Carlo simulation method for sampling the equilibrium ensemble of a lattice gas is obtained by 
considering a long series of small, random changes in the concentration that obey detailed balance. Before considering 
particular algorithms, it is necessary to factor the transition probability into the probability of attempting a move, 
A{c c'), and the probability of accepting the move V{c — > c'), i.e. TZ{c c') = A{c c')V{c c'). 

A commonly used and intuitively appealing Monte Carlo algorithm is the Metropolis algorithm. Given a particular 
configuration, all possible changes to the configuration are given the same attempt probability A. The attempted 
move is then accepted with probability 



■p(c c') = min exp 



Hic')-n{c) 



,1 , (14) 



RT 

which can be verified to obey the detailed-balance condition, Eq. ([T^), for all cases. A similar algorithm, known by 
the name Glauber also obeys detailed balance. It accepts the move with probability 

^^"^ ^'^ " 1 + exp [{n{c') ~ n{c)) /RT] ■ ^^^^ 

The two algorithms become identical as T — 0, except for H(c') — H(c) = 0. 

Another importance-sampling algorithm is useful when more than two configurations are possible for the portion 
of the surface being updated. For one move in this heat hath algorithm, the Boltzmann factor for every possible 
configuration is calculated, and these weights are used to randomly choose the final configuration. The new con- 
figuration is accepted automatically, with the initial configuration as one of the possible choices. In the case only 
two configurations are possible, the heat-bath algorithm reduces to the Glauber algorithm That this algorithm 
satisfies detailed balance can be easily verified. 

Commonly, Monte Carlo updates change only one adsorption site. However, in many applications equilibration can 
be slow, either because of the large fiuctuations that occur near second-order phase transitions (critical slowing-down 
|^,Q) or because of large local energy barriers. In such cases cluster update algorithms that involve several 
related sites in one Monte Carlo update can be useful. A specific example of a heat-bath cluster algorithm designed 



to overcome slowing-down due to high local barriers against single-particle moves is discussed in Sec. |VC. 

When importance sampling is used to incorporate the weighting, many thermodynamic averages can be estimated 
from the moments of the sampled quantities. For example, the g-th moment of the number of adsorbed particles of 
species X is 

s / N \ 9 

^^f-^-^E E^n ' (16) 
j=i \i=i / j 

where s is the number of samples taken during the simulation. The coverage is then just 

_ 1 OlnZ _ Ml 
^ Ndiflx/RT) iV ' ^ ' 

with the number of adsorption sites. The response functions, dQx/dp-x, which are related to the quasi-equilibrium 
CV current (see Sec. |l|), are p3[ | 



(18) 



djix RTd(fix/RTY 

s Mf - {MfY 



s~l N RT 

The response function is proportional to the variance of the number of adsorbed X particles, which corresponds to the 
fluctuations in the coverage. Equation (111) is an example of the class of important results in statistical mechanics. 
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called fluctuation-dissipation theorems |3j|. In our practical context it provides a convenient way to calculate the 
response functions without numerical differentiation. The cross- response function, dQx/dfiY^ is calculated entirely 
analogously, using the covariance of X and Y. 

In thermodynamic regions far from any phase transitions, the fluctuations in the numbers of adsorbate particles 
have a small correlation length. The surface can therefore be divided into independent pieces, and the fluctuations 
in adparticle numbers are proportional to the surface area. As a consequence the relative fluctuations, which are 
proportional to \/90x7^Ax/0x vanish with system size as N~^^'^. Near a second-order phase transition on the other 
hand, the correlation length becomes comparable to the system size, and the fluctuations and response functions 
depend on iV as power laws. The power-law exponents depend on the nature of the phase transition, and can be 
studied by finite-size scaling techniques |T^ , [l9| , |35| . 

The actual execution of a Monte Carlo simulation involves several steps. An initial configuration is constructed, 
followed by a sufficient number of Monte Carlo updates to allow it to reach equilibrium. A good habit is to construct 
the initial configuration using different phases for the right and left halves of the simulated system. This reduces the 
probability of relaxing to a long-lived metastable phase rather than true equilibrium. Also, several quantities should 
be sampled periodically during the equilibration process to ensure that the system has indeed reached thermodynamic 
equilibrium before equilibrium sampling begins. Equilibrium samples should be spaced far enough apart to be statis- 
tically independent. This is normally not a problem in parameter regions far from any phase transitions. Near phase 
transitions, however, critical slowing-down leads to very strong correlations between subsequent configurations and 
requires the sampling to be performed at longer intervals. While sampling with intervals that are too short leads one 
to underestimate fluctuations and thereby response functions, sampling too infrequently is a waste of computer time. 
Optimal sampling can be achieved by calculating time correlation functions |^ of the quantities being sampled and 
adjusting the interval accordingly. To improve accuracy, averages should also be performed over several trials, each 
consisting of equilibration followed by sampling. 

Finally, it must be noted that there is a difference between a single attempted move, often called a Monte Carlo 
step (mcs), and the more frequently encountered Monte Carlo step per site (MCSS). The later is useful because it 
is independent of the number of lattice sites and represents the smallest unit of updates over which the lattice as a 
whole can relax. These units should be distinguished clearly when reporting the results of Monte Carlo simulations. 



B. Dynamic Monte Carlo Simulations 



In the previous Subsection we discussed the basic interpretation of Monte Carlo simulation as a method to measure 
the equilibrium partition function and moments of quantities of interest. When it is used in this way, the simulation 
says nothing about the time dependence of the system. Consequently, one can use any algorithm that speeds up 
equilibration and reduces the correlations between subsequent configurations. The only requirement is that it must 
obey detailed balance in order to ensure that thermodynamic equilibrium is indeed reached and importance sampling 
can be used. 

A very different application of Monte Carlo simulation results when the number of MCSS is viewed as time. Then 
the simulation can be viewed as a "movie" of the dynamic process. Several examples of such movies can be found at 
the authors' Web site p6[ . 

It is clear that a particular Monte Carlo algorithm does not necessarily say anything about the time evolution of 
the corresponding chemical or physical system unless it is specifically constructed to do so. It is possible to construct 
a Monte Carlo simulation that evolves like the physical system, at least for dynamics that can be reasonably modeled 
as a stochastic process. This is the case when the dynamics are dominated by thermally activated processes such 
as adsorption, desorption, and diffusion. Here we describe the recipe for transforming a lattice-gas Hamiltonian, 
hopefully determined to reasonable accuracy by equilibrium studies, into such a dynamical model. 

Essentially, every mcs should correspond to a thermally activated transition between states separated by energy 
barriers, with the transition rates described by Arrhenius rates involving the barrier heights. These barrier heights 
can be constructed by a method essentially equivalent to the Butler- Volmer approximation |^^. This is done in the 
framework of activated-complex theory, or transition-state theory, with the initial and final states described by the 
lattice gas. The activated complex cannot be represented by the lattice gas, but a free energy can be assigned to it 
using the model we now describe. 

A schematic free-energy surface is presented in Fig. ^ for several possible lattice-gas configurations, '^c. First 
consider the two states with free energy 7i(^c) = Ti-i^c) = (depicted as solid line segments). The horizontal line 
segments correspond to the lattice-gas states, which are locally stable. The free energy is taken to increase linearly 
as the system moves toward the activated complex, which occurs where the two free-energy surfaces intersect. This 
linear assumption for the free-energy increase leads directly to the Butler- Volmer equation |27|. The height of the 
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free-energy barrier associated with this case, A, is a parameter which can be adjusted to describe particular processes. 

In general the slope of the free energy can also be adjusted to describe specific chemical systems. If the two slopes 
are described by the angles 9 and (/>, as shown in Fig. ||, then changes in the activated-complex free energy associated 
with 7i(c) ^ 0, can be found in terms of the symmetry factor 

ta-^ (19) 



tan 9 + tan (f) 

For a transition from the 'H(^c) = state to the HC^c) ^ state, the free energy of the activated complex is 

n* (^c, ^c) = A + (1 - a) H (^c) . (20) 

When a ^ 1/2, Fig. ^ is not symmetric about the activated complex. So if 7i('^c) ^ while H-i^c) = 0, the 
activated-complex free energy is 

n* fc,2c) = A + a7^ (3c) , (21) 

where the symmetry when a — 1/2 is obvious. For the most general case with TH^c) and 7i('*c), the free energy of 
the activated complex is 

n* (^c, ^c) = A + (1 - a) H (^c) + an (^c) . (22) 

This free energy is unaltered by reflection around the transition state. In addition, the energy barriers Ti* — Ti, (^c) 
and H* — Ti. (^c) are not altered by shifting the zero of the free-energy scale. 
The symmetric case, a = 1/2, gives the convenient relation 

n* (c, c') - A + i [t^ (c) + n (c') ] . (23) 

In this case the activated-complex free energy is found by adding a constant barrier to the arithmetic mean of the free 
energies of the lattice-gas states |Q. Once a consistent barrier scheme has been identified, a number of reasonable 
transition rates can be constructed, subject only to the requirement of detailed balance. The simplest rate is probably 
the one given by 

«,..,..„„,(_ -)»p(-?M_^), 

where vq is the attempt frequency, typically on the order of a phonon frequency (10^ - 10^^ Hz). Another useful 
choice is the Transition Dynamics Algorithm [Q, 

T^ic^c)^ ^ g(W.(c,c')--H(c))/i?T] g(W(c')--H-(c,c'))/fl,T] ' (^^) 

which is simply the product of the Glauber transition rates from c to the activated complex, and from the activated 
complex to c'. 

The single-ion processes included in the dynamical model discussed here are diffusion to nearest-neighbor sites and 
adsorption/desorption at single sites. We assume a single value for the angle in Fig. ^ for all configurations with the 
moving ion near the surface, which implies a = 1/2 for diffusion. In general, the angle associated with an ion far 
from the surface has a different value, but we make the reasonable ||2^ choice a — 1/2 so Eq. (^3|) can be used for 
all microscopic processes. On the other hand, we use a different barrier, A^, for each combination of ion and process 
P ( P = a for adsorption/desorption and P = d for diffusion). The values of the A^ determine the relative weights 
among the different microscopic processes. Reasonable agreement with experimental time scales would indicate that 
the combination of attempt frequencies and barrier heights chosen is realistic. 



C. The n-fold Way 



In many situations of experimental interest, a direct implementation of an importance-sampling Monte Carlo 
algorithm with Butler- Volmer transition rates would require impossibly long simulation times. The basic reason for 
this is that only a very small fraction of attempted Monte Carlo updates may actually result in a change of the current 
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configuration. Since we are now interested in the dynamical details, we cannot solve this problem by constructing a 
cluster algorithm, since this would generally change the dynamics completely. What is needed is an algorithm which 
can simply reduce the computer time needed to perform the large number of unsuccessful attempts that would be made 
between each successful update while remaining faithful to the original dynamics. A class of algorithms that perform 
this feat are known as rejection-free algorithms [ ^9|j40[ |. Many such algorithms have been invented and re- invented 
over the years (see Ref. |0|), they are all essentially equivalent. In these algorithms the acceptance probability is 
unity, but the probability of attempting a move is based on its Boltzmann weight. They can also be thought of as 
a Metropolis or heat-bath simulation with the number of rejected moves before the next successful move calculated 
analytically. The n-fold way introduced by Bortz, Kalos, and Lebowitz, is one of the best known and earliest of 
these methods. 

Rejection- free algorithms are implemented by keeping a list of all the possible moves. Each move is weighted by its 
rate TZ, so that making a random choice from the list is biased towards faster processes. The stochastic nature of real 
processes is captured by updating the simulation clock after each move by a random number that represents the time 
elapsed since the last move. For the current configuration, the overall rate for the system to have any move accepted 
is the sum of all the individual move rates, TZc- Since the probability of escaping from the current configuration has 
a constant rate, the elapsed time until the accepted change has an exponential distribution, and can be found using 

At = -7^^l In (r) , (26) 

where r is a random number uniformly distributed on (0, 1). Here At is a continuous time in units of mcs. 

The weighted list can be visualized as a line segment whose length equals the sum of all the weights. Then, choosing 
from the weighted list is accomplished by picking a random point along the line segment, i.e. picking a number on 
the interval [0,1) and multiplying by the segment length. Determining which choice this corresponds to involves 
searching the list of process rates. This searching can become an important bottle neck in the implementation of 
n-fold way algorithms, but the effect can be reduced by using a binary tree to search for the choice from the weighted 
list. Building a binary tree of partial sums (of the rates) to speed up the search is not difficult, and a great deal of 
help is available from many introductory texts on data structures [Q. The rates make up the leaves at the bottom 
of the tree, the value of each node is the sum of all the leaves below it. As the simulation proceeds, the values of the 
leaves change and the nodes must be updated. Two approaches can be taken. The simplest is to move up the tree 
from the changed leaf, explicitly adding the values of the children to find the value of the node. A faster approach 
is to calculate the change in the leaf value, and add this to all the nodes above the leaf Unfortunately, with this 
faster approach truncation errors will eventually cause the nodes to not accurately reflect the values of the leaves. 
Practically this method can be used, but the entire tree should be recalculated every hundred thousand MCSS or so. 



IV. SPECIFIC EXAMPLES 



A. Br on Ag(lOO) 



Experimental CVs of bromine adsorption onto silver (100) single-crystal surfaces display a broad peak followed by a 
stronger, sharper peak at more positive potentials [n3|42]. In situ surface x-ray scattering shows that the sharp peak 
corresponds to a second-order, i.e. continuous, phase transition from a disordered low-coverage phase at negative 
potentials to an ordered c(2 x 2) phase at more positive potentials jl^]. Since the adsorption sites for the bromine are 
the four- fold hollow sites of the square Ag(lOO) surface |^,^, it is natural to model this system with a lattice-gas 
Hamiltonian of the form given by Eq. (|l|), using a square lattice of adsorption sites. 

The ordered phase corresponds to a maximum bromine coverage of O = 0.5. This is related to the adsorbed bromine 
being bigger than the Ag(lOO) lattice spacing, such that adsorption in nearest-neighbor hollow sites is very unfavorable 
energetically. This packing feature can be modeled as an infinitely strong repulsion between nearest-neighbor bromine, 
<I>(^) = — oo. Relaxing the nearest-neighbor repulsion to finite, but large, values does not significantly improve the 
agreement with experiment [Q. For separations larger than this, pair interaction energies are modeled by unscreened 
dipole-dipole interactions: 

= — , (27) 

r[n) 

where r{n) is the distance between sites, measured in units of the lattice constant. This interaction falls off reasonably 
rapidly with distance, and only interactions for n < 13 were considered in the simulations discussed here. This 
corresponds to approximately 87% of the total dipole-interaction energy. Other lattice-gas simulations llj] used 
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screened dipole-dipole interactions, but the effect of screening is relatively small for Br~ at lengths less than the 
truncation used here. 

If (I)(2) = 0, the only interaction is the infinite nearest-neighbor repulsion, and the lattice-gas model is just a hard- 
square model G^jH. The hard-square model has a critical point at /I « 1.334i?T, which corresponds to the onset 
of percolation |44| at a critical coverage 8c ~ 0.37 Below this concentration a disordered gas of adsorbate 

particles forms, but above Qc an ordered c(2 x 2) phase with coverage near 1/2 forms. Typical configurations of the 
two phases are illustrated in Fig. |^. Numerical studies of the hard-square model with repulsive next-nearest neighbor 
interactions indicate that a 8 = 1/4 phase (which is not seen experimentally in the system studied here) also 

becomes important for strongly repulsive This constrains the repulsive interactions in the lattice-gas model. If 

the interactions are too strong, a plateau in the isotherm occurs around 8 = 1/4, in disagreement with experiments. 

To generate a numerical adsorption isotherm, the bromine coverage at equilibrium was measured from Monte Carlo 
simulations at a discrete set of electrochemical potentials using a simple Metropolis algorithm. Since no finite-size 
scaling analysis was planned, only one lattice size was used: a 32x32 system with periodic boundary conditions. The 
size and boundary conditions used illustrate general rules for simulation of lattice-gas models: 

• The system size must be chosen as a multiple of the periodicities of all important ordered phases. In the present 
case this is easy, since the c(2x2) phase has period two in both directions. 

• Unless there is a specific reason to do otherwise, one eliminates edge and corner effects by using periodic 
boundary conditions. Geometrically, this can be visualized as wrapping the system onto a torus. The system 
size should be larger than the longest-ranged interaction to avoid wrapping all the way around. 



• Near critical points, fluctuations are correlated over lengths comparable to the system size (see Sec. III). Accurate 
simulation near critical points therefore requires large systems. 

Some sophistication is needed in the implementation to make the simulation reasonably efficient. The long-range 
repulsion is truncated at the 13-th neighbor, a separation of 5 lattice constants. Due to the long interaction range, a 
significant amount of computer time is spent evaluating the energy change due to an attempted move. To this end, a 
data structure tracking the number of n-th neighbor bromines for each site is maintained and updated only after each 
successful move. For this simulation the Metropolis acceptance rate is around 10%, so this leads to a considerable 
speed-up. A more efficient simulation using a heat-bath algorithm with cluster updates is currently being investigated. 

The experimental isotherm, determined by chronocoulometry p3| , p^ for 0.1 mM Br, is shown in Fig. | (circles), 
along with two theoretical results. The effective electrovalence, z, can be determined from the experimental isotherms 
by assuming the coverage depends only on the electrochemical potential and then fitting z to isotherms for different 
electrolyte concentrations using Eq. (g). Using the experimental results from Ref. |l^ we find z « —0.736, which is 
in good agreement with other estimates [Q^. With z fixed at this value and <I>*^^-' taken as a fitting parameter, the 
best fits for both the single-sublattice Frumkin isotherm |l5j (dotted line) and the full lattice-gas model evaluated by 
Monte Carlo simulations give ^'^^^ « — 24meV « — 2.3kJ/mol. 

Both theoretical isotherms agree qualitatively with the experimental data, but there is a noticeable difference near 
8 « 0.35 where the lattice-gas model shows the singularity in the slope of the isotherm associated with the phase 
transition. The Frumkin isotherm does not contain this singularity, but only a gradual cross-over from low to high 
coverage. In this region, there is considerably better agreement between the experimental and Monte Carlo results. 
Reproducing this behavior is important, because the singularity gives rise to a peak in the quasi-equilibrium CV 
(not shown here) which corresponds to the peak seen in experiments. Thus, the Monte Carlo simulations provide a 
theoretical description of the essential part of this electrochemical adsorption process that the mean-field Frumkin 
isotherm cannot. 

B. Cu on Au(lll) in the Presence of Sulfate 

In underpotential deposition (UPD), a monolayer or less of one metal is electrochemically adsorbed onto another 
in a range of electrode potentials more positive than those where bulk deposition occurs . The UPD of copper on 
Au(lll) electrodes in sulfate-containing electrolytes has been intensively studied, both experimentally and theoreti- 
cally. A discussion of the literature until 1995 is included in Ref. |]l9|. The most striking feature in CV experiments on 
this system is the appearance of two peaks, separated by ~ 50-150 mV, upon addition of Cu^+ ions to the sulfuric-acid 
electrolyte [Esl-lsTl. Typical simulation CV profiles, estimated by both quasi-equilibrium and dynamic methods, are 



shown in Fig. |^ In the potential range between the peaks, the adsorbate layer is believed to have a (v^ x v^) 
structure consisting of 2/3 monolayer (ML) copper and 1/3 ML sulfate |lj]. This structure, which is illustrated in 
Fig. [|, is strongly supported by in situ x-ray diffraction data |5^ . 
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The two-component lattice-gas model used in references 53|, is a refinement of the original Huckaby and Blum 
model for this system. It was extended with interactions of intermediate range, similar to the extensions to the 
original model discussed in Ref. [Q. These models are based on the assumption that the sulfate coordinates the 
triangular Au(lll) surface through three of its oxygen atoms, with the fourth S-0 bond pointing away from the 
surface. This gives the sulfate a triangular footprint with a 0-0 distance that reasonably matches the lattice constant 
of the triangular Au(lll) unit cell. The copper is assumed to compete for the same adsorption sites as the sulfate. 
The configuration energies are given by Eq. (^) with A — C for copper and B = S for sulfate. 

Just as for the Br/Ag(100) system discussed in Sec. IV A , the sulfate ion is sufficiently large that steric interactions 
prevent two sulfates from binding simultaneously to nearest-neighbor sites. This is modeled as an infinite repulsion 
between nearest-neighbor sulfates, $33* — —00. Neglecting other interactions, this corresponds to the hard-hexagon 
model ||l7| , |l8| , ^ , ^ which can be solved analytically [^,0. In the limit of strong sulfate adsorption, //g +00, 
the infinite nearest-neighbor repulsion yields a 1/3 ML sulfate phase. Denoting a phase with a, XxY unit cell as 

(XxY)q^, the highest density sulfate-only phase is (VSx-v/S)^^- 

The (rather complicated) complete ground-state diagram for a lattice gas model of copper UPD is shown in Fig. [|. 

Experimentally, the two CV peaks separate a copper monolayer, an ordered (\/3 x ^8)^/3 phase, and a disordered 

low-coverage (Ix 1)q phase. In experiments involving sulfate but no copper, the (1x1)q and (V3xV7)o^^ phases are 
observed ||5^,Q. The desire to reproduce this sulfatc-only phase largely dictated the values of the effective sulfate- 
sulfate interactions used in the model ||l^. In the model, more positive electrode potentials cause the sulfate to form 
its saturated {VSxy/S)^^^ hard-hexagon phase. However, in reality other chemical processes, such as surface oxidation, 
probably occur before this phase is fully achieved. 

The effective electrovalences, zg and zc, must be determined from experiments [^7|-^. Zhang, et ai, obtained 
the values, zq ~ -fl.7 and zs ~ —1-1, using data for the dependence of the CV-peak separation on electrolyte 
concentration from Omar et al. [^7|| . These values are consistent with independent estimates [ 58| , |59| . As explained in 
Sec. g constant-concentration isotherms are lines of slope zq/zc in the phase diagram. The isotherm used for this 
discussion is given by the dashed fine in Fig. for 100 < £; < 300 mV vs Ag/AgCl. 

Starting from the negative end, we scan in the direction of positive electrode potential (upper left to lower right 
in the figure). Near the CV peak, at approximately 185 mV vs Ag/AgCl, the sulfate begins to compete with copper 
for the gold surface sites, resulting in a third of the copper desorbing into the bulk and being replaced by sulfate. 
Due to the strong effective attraction between the copper and sulfate adparticles, a mixed (-s/Sx V^)^^^ phase is 
formed (see Fig. |]), which extends through the entire potential region between the two CV peaks. As the CV peak 
at approximately 250 mV is reached, most of the copper is desorbed within a narrow potential range. As it is thus 
deprived of the stabilizing infiuence of the coadsorbed copper, the sulfate is partly desorbed, reducing Os from 1/3 
to approximately 0.05. 

Both equilibrium and time-dependent Monte Carlo simulations have been performed for this system, using the 
same model parameters as in Ref. |l^. The equilibrium simulations, which produced the quasi-equilibrium CV 
currents shown in Fig. |^, were performed with a heat-bath algorithm. To circumvent energy barriers that reduce the 
acceptance rate for single-site updates, a cluster-update algorithm treating two nearest-neighbor sites simultaneously 
was employed. The simulations were performed on a rhombus-shaped triangular-lattice system of size 30 x 30 with 
periodic boundary conditions. At each value of the electrode potential, the system was equilibrated for between 50 
and 500 MCSS, after which sampling was performed every 5 MCSS for up to 5000 MCSS. Averaging over 20 trials 
was performed, except near the sharp peak, where 215 trials were conducted. 

Dynamic simulations were also performed on a 30 x 30 system, using an n-fold way algorithm which includes 
adsorption, desorption, and nearest-neighbor diffusion moves for both Cu and sulfate. We chose = 10^° s~^ in 
Eq. ( p^ ) and used the energy barriers = = 35kJ/mol, Ag = 55kJ/mol, and A^ = 50kJ/mol. The resuh for 
a dynamic simulation of a CV with scan rate dE/dt = lOmV/s are also included in Fig. ||. After initial equilibration 
at a constant electrode potential and room temperature for 1 sec, the dynamic simulation proceeds by changing E in 
steps of 1 mV, simulating for a time determined by the scan rate, and then sampling once. This process is continued 
until the end of the scan is reached. The results from each potential scan are averaged over 75 independent simulation 
runs, and separate simulations are used for the positive-going and negative-going scans. The dynamic CV has broader 
peaks than the quasi-equilibrium CV obtained from the equilibrium simulation. The peaks are also displaced in the 
scan direction, relative to the quasi-equilibrium peak positions. Clearly, dynamic simulations provide information 
about possible dynamic effects in CV experiments with finite potential-scan rates. 

Dynamic Monte Carlo simulations can also be used to study the current transients that occur after the electrode 
potential is rapidly stepped across a transition | ]60[ |. In this case, the system is equilibrated on one side of the phase 
transition (located by the peak positions in the quasi-equilibrium CV), then the potential is changed instantaneously 

across the transition to its final value. The results for potential steps simulations across the (VSx \/3)2/3 ^ (Ix 1)q 
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transition are shown in Fig. |g. 

Fig. 6a records current transients for potential steps that start in the ordered phase and end between 24 and 33 mV 
past the transition on the low-coverage side. The first current transient corresponds to the relatively fast relaxation 
to a metastable state with slightly less adsorbed Cu. This is followed be a second current maximum as domains of 
empty surface nucleate and grow. The trend of the second maximum to early times and larger currents as the step 
size is increased agrees with recent experiments |l60| . 

The current transient for negative steps in the cell potential is shown in Fig. 6b. The transient simply decreases with 
time, corresponding to adsorption onto the empty surface that does not require nucleation, i.e. the surface is unstable 
with respect to adsorption after the potential step. This difference between the observed transients after positive- 
and negative-going potential steps across this particular transition is precisely the difference seen in experiments of 
Cu with SO4 adsorbing onto Au(lll) |60|. 



C. Urea on Pt(lOO) 

The adsorption of urea onto Pt(lOO) surfaces requires an even more complicated lattice-gas model than the two 
systems discussed above. Experimental and simulated quasi-equilibrium CV currents for this system are shown 
together in Fig. 0. The urea coverage Qtj, measured in situ P,|2^, changes over a narrow potential range from near 
zero on the negative side of the CV peak to approximately 1/4 ML on the positive side. Ex situ LEED at potentials 
on the positive side of the peak show an ordered c(2 x 4) structure, consistent with a urea coverage of 1/4 ML. On 
the negative side of the peak, only an unreconstructed (1x1) surface with 8h = 1 is seen (H denotes hydrogen). 

The model developed to account for these observations is based on the assumption that urea, CO(NH2)2, coordinates 
the platinum through its nitrogen atoms, with the C=0 group pointing away from the surface. Since the unstrained 
N-N distance in urea matches the lattice constant of the square Pt(lOO) surface quite well, urea is assumed to occupy 
two nearest-neighbor adsorption sites on the square Pt(lOO) lattice. Coulometry data indicate that the hydrogen 
coverage changes from near unity on the negative-potential side of the CV peak to near zero over the same potential 
range where urea becomes adsorbed. Therefore, hydrogen is assumed to adsorb in the same on-top positions as 
the urea nitrogen atoms. In the resulting model hydrogen is adsorbed at the nodes and urea on the bonds of a 
square lattice. Simultaneous occupation by two or more urea molecules of bonds that share a node is excluded, as is 
occupation by hydrogen of a node adjacent to a bond occupied by urea. An illustrative Monte Carlo snapshot of a 
typical c(2 x 4) configuration with thermal fluctuations is shown in Fig. 0. 

The configuration energies of this lattice-gas model are given by Eq. ^) with A=U (urea) and B=H (hydrogen). 
The effective lattice-gas interactions were determined from ground-state calculations followed by numerical Monte 
Carlo simulations of quasi-equilibrium CV currents ||6|,p2[|. In order to stabilize the observed c(2 x 4) phase, effective 
interactions were included through n = 8 P,p2[. 

The numerical Monte Carlo simulations, which used systems with up to 32 x 32 square-lattice unit cells, were 
performed with a heat-bath cluster algorithm. The clusters are cross-shaped regions, consisting of five nearest-neighbor 
nodes plus their four connecting bonds. After symmetry reductions, these clusters have 64 different configurations. As 
a result, the code is relatively slow, measured in CPU time per attempted update. However, the algorithm effectively 
removes energy barriers due to "jamming" configurations, especially deep in the c(2 x 2) phase. In that phase region 
the cross-cluster algorithm outperforms an algorithm based on a minimal cluster consisting of one bond plus its 
two adjacent nodes by at least a factor of several million in terms of CPU time to achieve equilibration. The good 
agreement of the quasi-equilibrium CV for the lattice-gas model with the experimental results, which is shown in 
Fig. again demonstrates the utility of Monte Carlo simulations in the study of electrochemical adsorption. 



V. SUMMARY 



Our goal has been to convince practitioners of interfacial electrochemistry that Monte Carlo simulation is a viable 
and useful research tool, well within the computational means of most laboratories. We have described how numerical 
Monte Carlo simulations of lattice-gas models for specific electrochemical adsorption can provide realistic estimates 
for a number of experimentally observable quantities. These include strictly equilibrium ones such as coverages and 
surface charges, quasi-equilibrium quantities like CV currents for extremely slow potential-scan rates, and genuine 
nonequilibrium quantities such as current transients and CV currents at high scan rates. Information about model 
parameters, including effective lateral interaction energies and diffusion barriers for adsorption/desorption and lateral 
diffusion, can also be obtained. In this Chapter we have presented detailed discussions of three specific systems. In 
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doing so, we have attempted to give enough detail about the simulations and references to general texts on simulation 
that it should enable the reader to write and use her or his own Monte Carlo simulations. 

Most models with realistic geometric structure and interactions are not exactly solvable, and the solution of even 
highly simplified models must be accomplished by approximate means. Monte Carlo simulation represents a superior 
alternative to most mean-field based approaches because it provides a much more accurate connection between the 
model parameters and the values of the calculated observables. If comparison of experimental and theoretical results is 
going to yield reliable information about the values of microscopic model parameters, Monte Carlo simulation should 
usually be the method of choice. 
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FIG. 1. Ground-state diagram for the lattice-gas model of copper UPD on Au(lll), shown in the (/Is, Ac) plane. The 
subscripts S and C denote sulfate and Cu, respectively. The effective interactions are given in Ref. flSl]. The solid lines 
represent zero-temperature phase boundaries, the dotted line is a zero-temperature phase boundary that is not observed in 
room-temperature simulations. The voltammetric scan path for the simulations discussed here is represented as the dashed 
line. Its end points correspond to £'=110 mV (upper left) and 300 mV (lower right) vs Ag/AgCl, respectively. The phases are 
denoted as (XxY)®^. The ground-state configurations corresponding to the main phases are also shown. The adsorption sites 
are shown as o, and Cu and sulfate are denoted by • and A, respectively. (Adapted from Ref. VM.) 



FIG. 2. Schematic free-energy surfaces for thermally activated transitions in the dynamic Monte Carlo model. The text 
describes how these free-energy surfaces, based on the Butler- Volmer model, are used to calculate energy barriers for Monte 
Carlo processes. In our simulations processes corresponding to adsorption, desorption, and surface diffusion are included. 



FIG. 3. Snapshots of configurations from simulation of bromine adsorption on Ag(lOO) for (a) the low-coverage disordered 
phase and (b) the ordered c(2 x 2) phase. The light gray circles represent surface Ag atoms, and the black circles the adsorbed 
Br atoms. The Br adsorb to the four-fold hollow sites, which form a square lattice. 
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FIG. 4. Experimental and theoretical isotherms for bromine adsorption on Ag(lOO). The circles are chronocoulometry results 
from Ref. The solid curve is the numerical Monte Carlo result, while the dashed curve is the best-fit single-sublattice 

Frumkin isotherm. The Monte Carlo isotherm displays the singularity associated with the phase transition, while the Frumkin 
isotherm does not. 



FIG. 5. Simulated CV profiles for copper UPD on Au(lll), corresponding to the scan path in Figure Solid curves are 
dynamic simulations with scan rate lOmV/s, while dotted curves are quasi-equilibrium estimates. The model parameters are 
the same as in Ref. |l9|l. 



FIG. 6. Simulated current transients after potential steps across the (Ix 1)q ^ (\/3x^/3)2/3 transition, (a) Positive steps, 
(b) Negative step. The observed asymmetry is the same in the simulations as in experiments. (Reproduced with permission 
from Ref. [^. Copyright 1998 Elsevier Science.) 



FIG. 7. Room-temperature CV profiles for urea on Pt(lOO) in 0.1 M HCIO4. Experimental (dashed curves) and simulated 
(O and solid curve) quasi-equilibrium normalized CV currents, i/{dE/dt), in elementary charges per mV per Pt(lOO) unit cell, 
at 1.0 mM bulk urea. The dashed curves are representative negative-going voltammograms; two are shown to demonstrate 
variations between individual measurements. The effective interactions and electrovalences are given in the caption of Fig. 4 of 



Ref. 
Ref. 



2|. The simulation results shown here are for a system of 32 x 32 Pt(lOO) unit cells. (Reproduced with permission from 
|. Copyright 1995 Elsevier Science.) 



FIG. 8. Left; The microscopic model for urea on Pt(lOO), showing urea molecules with their NH2 groups (dark gray) on the 
square Pt(lOO) lattice sites and their C=0 groups (C black and O lighter gray) pointing away from the surface. The hydrogen 
atoms are shown as single light gray spheres. Right: A typical equilibrium configuration generated by Monte Carlo in the 
ordered c(2 x 4) phase at — 60mV, using the model parameters given in Ref. [^. Urea molecules are shown as filled rectangles 
on the lattice bonds and hydrogen as • at the nodes. (Reproduced with permission from Ref. Q. Copyright 1995 Elsevier 
Science.) 
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